#updated 1/8/2015 to work with R 3.1.2 and lme4 1.1-7
#==================================================================
#Tall Thistle: data not available for Seed production
#data available from Leland see below
#==================================================================



#==================================================================
#Bull Thistle: 
#first calculate probability of seeds >0
#then calculate seed production for plants that produce >0 seed
#==================================================================


BTseeds=read.csv("Reproduction 2006&2007.csv",header = T,na.strings=c("NA","na"))
#converting variables in "bolt" data to numeric/categorical variables
BTseeds$comp = as.factor(BTseeds$comp)
BTseeds$herb = as.factor(BTseeds$herb)
BTseeds$site = as.factor(BTseeds$site)
BTseeds$year = as.factor(BTseeds$year)
BTseeds$site.year = paste(BTseeds$site,BTseeds$year,sep=".")
#BTseeds$noseeds=round(BTseeds$noseeds)
hist(BTseeds$noseeds)
table(BTseeds$comp,BTseeds$herb)

library(lme4)

# relevel competition factor
#BTseeds$comp = factor(BTseeds$comp,levels=c("Low","Medium","High"))
BTseeds$herb = factor(BTseeds$herb,labels=c("Reduced","Ambient"))
BTseeds = BTseeds[complete.cases(BTseeds[,14]),]
plot(noseeds~biomass, data = BTseeds[BTseeds$year==2006,])
points(noseeds~biomass, data = BTseeds[BTseeds$year==2007,],col = "red")
#BTseeds$scaleseeds = scale(BTseeds$noseeds)
BTseeds$logseeds = log(BTseeds$noseeds+1)
BTseeds$logbiomass = log(BTseeds$biomass)

table(BTseeds$comp[BTseeds$noseeds==0],BTseeds$herb[BTseeds$noseeds==0]) #zero seed production only for ambient herbivory
#Probability of zero seed production under reduced herbivory is zero
#So I calculate this probability only for ambient herbivory 
AmbientSeeds=BTseeds[BTseeds$herb=="Ambient",]

Ambient.mm0 = glmer(seedproducedyesno~logbiomass + comp +  
		logbiomass:comp + 
		(1|site.year),  family=binomial,data=AmbientSeeds)
Ambient.mm1 = glmer(seedproducedyesno~logbiomass + comp + 
		logbiomass:comp + 
		(comp|site.year),  family=binomial,data=AmbientSeeds)
sapply(list(Ambient.mm0,Ambient.mm1),BIC)
#BIC is smallest in Ambient.mm0, so I use (1|site.year)

AmbientSeeds.models = list("~logbiomass + comp + logbiomass:comp",
		"~logbiomass + comp ",
		"~comp",
		"~logbiomass",
		"~1")
AmbientSeeds.fits = lapply(AmbientSeeds.models,function(ff)glmer(paste("seedproducedyesno",ff,"+ (1|site.year)"),data=AmbientSeeds,family=binomial))
bic = unlist(sapply(AmbientSeeds.fits,BIC))
dbic = bic-min(bic)
wbic = exp(-dbic/2)/sum(exp(-dbic/2))
data.frame(unlist(AmbientSeeds.models),dbic,wbic)[order(dbic),]
loglik=unlist(sapply(AmbientSeeds.fits,BIC))
#k=sapply(tg.fits,function(mm)length(fixef(mm)))+1
k=sapply(AmbientSeeds.fits,function(mm)length(fixef(mm)))+1

stat.res=data.frame(unlist(AmbientSeeds.models),k,loglik,dbic,wbic)[order(dbic),]
write.csv(stat.res,"ZeroSeeds_res_BT.csv")

#model 4 has 89% of the weight

# make a plot
nd = data.frame(logbiomass=rep(seq(-0.7,6.1,0.1),times=2),
		comp=factor(rep(c("0","1"),each=69),levels=levels(AmbientSeeds$comp)))

results = matrix(NA,nrow=nrow(nd),ncol=length(AmbientSeeds.fits))
for (i in 1:length(AmbientSeeds.fits)){
	# generate predictions from a glmer model using the data from inside the model, mostly. 
	fm = AmbientSeeds.fits[[i]]
	results[,i] = predict(fm,newdata=nd,re.form=~0)
}
# model averaged prediction for logbiomass parameter
nd[,"nonzero"] = apply(results,1,function(x)sum(x*wbic))

xyplot(nonzero~logbiomass|comp,data=nd,type="l",
panel=function(...){
		panel.xyplot(...)
		pn = which.packet()
		#ltext(0,0.5,pn[1])
		pick = AmbientSeeds$comp == levels(AmbientSeeds$comp)[pn[1]] 			
		ticks = AmbientSeeds[pick,c("logbiomass","seedproducedyesno")]
		panel.rug(ticks[ticks$seedproducedyesno==0,1])
		panel.rug(ticks[ticks$seedproducedyesno==1,1],y=NULL, regular=FALSE)
	},
ylim=c(0,1))

#------------------
coef.names = names(fixef(AmbientSeeds.fits[[1]]))
coef.results = matrix(NA,nrow=length(AmbientSeeds.fits),ncol=length(coef.names))
coef.se = matrix(NA,nrow=length(AmbientSeeds.fits),ncol=length(coef.names))
for (i in 1:length(AmbientSeeds.fits)){
	coef.results[i,] = sapply(coef.names,function(cf)fixef(AmbientSeeds.fits[[i]])[cf],USE.NAMES=FALSE)
	se = sqrt(diag(vcov(AmbientSeeds.fits[[i]])))
	names(se) = names(fixef(AmbientSeeds.fits[[i]]))
	coef.se[i,] = sapply(coef.names,function(cf)se[cf])
}

coef.avg = apply(coef.results,2,function(x)sum(x*wbic,na.rm=TRUE))
t1 = sweep(coef.results,2,coef.avg)
t2 = sqrt(coef.se^2 + t1^2)
coef.avg.se = apply(t2,2,function(x)sum(wbic*x,na.rm=TRUE))

nd.intercept = data.frame(logbiomass=rep(0,times=2),
		comp=factor(c("0","1")))
mm.intercept = model.matrix(as.formula(AmbientSeeds.models[[1]]),data=nd.intercept)

nd.slope = data.frame(logbiomass=rep(1,times=2),
		comp=factor(c("0","1")))

mm.slope = model.matrix(as.formula(AmbientSeeds.models[[1]]),data=nd.slope)

mm.slope[,c(1,3)] = 0

int.results = array(NA,dim=c(length(AmbientSeeds.fits),nrow(mm.intercept)))
int.se = matrix(NA,nrow=length(AmbientSeeds.fits),ncol=nrow(mm.intercept))
slope.results = array(NA,dim=c(length(AmbientSeeds.fits),nrow(mm.slope)))
slope.se = matrix(NA,nrow=length(AmbientSeeds.fits),ncol=nrow(mm.intercept))

for (i in 1:length(AmbientSeeds.fits)){
	fef = rep(0,times=length(coef.names))
	names(fef) = coef.names
	vcf = matrix(0,nrow=length(coef.names),ncol=length(coef.names))
	rownames(vcf) = coef.names
	colnames(vcf) = coef.names

	fe = fixef(AmbientSeeds.fits[[i]])
	fef[names(fe)] = fe
	vcf[names(fe),names(fe)] = as.matrix(vcov(AmbientSeeds.fits[[i]]))
	int.results[i,] = mm.intercept %*% fef
	slope.results[i,] = mm.slope %*% fef
	for (j in 1:nrow(mm.intercept)){
		pick = which(mm.intercept[j,]>0)
		int.se[i,j] = sqrt(sum(vcf[pick,pick]))
		pick = which(mm.slope[j,]>0)
		slope.se[i,j] = sqrt(sum(vcf[pick,pick])) # make it the se not the variance
	}
	
}

int.avg = apply(int.results,2,function(x)sum(x*wbic,na.rm=TRUE))
t1 = sweep(int.results,2,int.avg)
t2 = sqrt(int.se^2 + t1^2)
int.avg.se = apply(t2,2,function(x)sum(wbic*x,na.rm=TRUE))

slope.avg = apply(slope.results,2,function(x)sum(x*wbic,na.rm=TRUE))
t1 = sweep(slope.results,2,slope.avg)
t2 = sqrt(slope.se^2 + t1^2)
slope.avg.se = apply(t2,2,function(x)sum(wbic*x,na.rm=TRUE))

#Probability of zero seed production under reduced herbivory is zero
#Thus I don't have a slope and intercept for reduced herbivory
IPM.parameters[5:6,23:26]=NA

#So I calculate this probability only for ambient herbivory 
#com = 0 refers to "Medium", and com = 1 refers to "High"
for (i in 2:3)  IPM.parameters$zero.seed.intercept[i]=int.avg[i-1]
for (i in 2:3)  IPM.parameters$zero.seed.intercept.se[i]=int.avg.se[i-1]
for (i in 2:3)  IPM.parameters$zero.seed.slope[i]=slope.avg[i-1]
for (i in 2:3)  IPM.parameters$zero.seed.slope.se[i]=slope.avg.se[i-1]




#===============================
#Remove zero seed production
library(lme4)
remove=which(BTseeds$noseeds==0)

#uses lmer if family = gausian
BTseeds.mm0=lmer(logseeds~logbiomass + comp + herb + 
		logbiomass:comp + logbiomass:herb + comp:herb + 
		(1|site.year),  data=BTseeds[-remove,])
BTseeds.mm1 = lmer(logseeds~logbiomass + comp + herb +
		logbiomass:comp + logbiomass:herb + comp:herb + 
		(herb|site.year) , data=BTseeds[-remove,])
BTseeds.mm2 = lmer(logseeds~logbiomass + comp + herb +
		logbiomass:comp + logbiomass:herb + comp:herb + 
		(comp|site.year),  data=BTseeds[-remove,])
# OUCH, we actually used AIC here, not BIC for selection of random effects
sapply(list(BTseeds.mm0,BTseeds.mm1,BTseeds.mm2),AIC)
Res1=resid(BTseeds.mm1)
plot(fitted(BTseeds.mm1),Res1)
#which(Res1< -40000)
#which(fitted(BTseeds.mm5) >80000) 
#Remove=c(56,142,172)
qqnorm(Res1)
qqline(Res1)

BTseeds.models = list("~logbiomass + herb + comp + logbiomass:comp + logbiomass:herb + herb:comp",
		"~logbiomass + herb + comp + logbiomass:comp + logbiomass:herb",
		"~logbiomass + herb + comp + logbiomass:comp + herb:comp",
		"~logbiomass + herb + comp + logbiomass:herb + herb:comp",
		"~logbiomass + herb + comp + logbiomass:herb",
		"~logbiomass + herb + comp + herb:comp",
		"~logbiomass + herb + comp + logbiomass:comp",
		"~logbiomass + herb + comp",
		"~logbiomass + herb +  logbiomass:herb",
		"~logbiomass + herb",
		"~logbiomass + comp + logbiomass:comp",
		"~logbiomass + comp",
		"~herb + comp + herb:comp",
		"~herb + comp",
		"~herb",
		"~comp",
		"~logbiomass",
		"~1")
BTseeds.fits = lapply(BTseeds.models,function(ff)lmer(paste("logseeds",ff,"+ (herb|site.year)"),data=BTseeds[-remove,]))
#BTseeds.fits = lapply(BTseeds.models,function(ff)glmer(paste("logseeds",ff,"+ (herb|site.year)"),data=BTseeds,family=gaussian))

bic = unlist(sapply(BTseeds.fits,BIC))
dbic = bic-min(bic)
wbic = exp(-dbic/2)/sum(exp(-dbic/2))
data.frame(unlist(BTseeds.models),dbic,wbic)[order(dbic),]
# model 10 has 35% of the weight, and model 9 has 32%, then come models 9 25%; those three modes have 96% of teh weight
loglik=unlist(sapply(BTseeds.fits,BIC))
k=sapply(BTseeds.fits,function(mm)length(fixef(mm)))+3
stat.res=data.frame(unlist(BTseeds.models),k,loglik,dbic,wbic)[order(dbic),]
write.csv(stat.res,"Seeda_res_BT.csv")

BTseeds.fits[[10]]

# make a plot
nd = data.frame(logbiomass=rep(seq(-0.7,6.7,0.1),times=4),
		comp=factor(rep(c(0,1),each=150),levels=levels(BTseeds$comp)),
		herb=factor(rep(rep(c("Reduced","Ambient"),2),each=75),levels=levels(BTseeds$herb)))

results = matrix(NA,nrow=nrow(nd),ncol=length(BTseeds.fits))
for (i in 1:length(BTseeds.fits)){
	# generate predictions from a glmer model using the data from inside the model, mostly. 
	fm = BTseeds.fits[[i]]
	results[,i] = predict(fm,newdata=nd,re.form=~0)
}
# model averaged prediction for veg.bio07 parameter
nd[,"noseeds"] = apply(results,1,function(x)sum(x*wbic))


xyplot(noseeds~logbiomass|comp+herb,data=nd,type="l",
panel=function(...){
		panel.xyplot(...)
		pn = which.packet()
		ltext(0,0.5,paste(pn,collapse=", "))
		pick = BTseeds$comp == levels(BTseeds$comp)[pn[1]] & 
			 BTseeds$herb == levels(BTseeds$herb)[pn[2]]
		ticks = BTseeds[pick,c("logbiomass","logseeds")]
		lpoints(ticks$logbiomass,ticks$logseeds)
	}
)


#------------------
coef.names = names(fixef(BTseeds.fits[[1]]))
coef.results = matrix(NA,nrow=length(BTseeds.fits),ncol=length(coef.names))
coef.se = matrix(NA,nrow=length(BTseeds.fits),ncol=length(coef.names))
for (i in 1:length(BTseeds.fits)){
	coef.results[i,] = sapply(coef.names,function(cf)fixef(BTseeds.fits[[i]])[cf],USE.NAMES=FALSE)
	se = sqrt(diag(vcov(BTseeds.fits[[i]])))
	names(se) = names(fixef(BTseeds.fits[[i]]))
	coef.se[i,] = sapply(coef.names,function(cf)se[cf])
}

coef.avg = apply(coef.results,2,function(x)sum(x*wbic,na.rm=TRUE))
t1 = sweep(coef.results,2,coef.avg)
t2 = sqrt(coef.se^2 + t1^2)
coef.avg.se = apply(t2,2,function(x)sum(wbic*x,na.rm=TRUE))

nd.intercept = data.frame(logbiomass=rep(0,times=4),
		comp=factor(rep(0:1,each=2)),herb=factor(rep(c("Ambient","Reduced"),2),levels=levels(BTseeds$herb)))

mm.intercept = model.matrix(as.formula(BTseeds.models[[1]]),data=nd.intercept)

nd.slope = data.frame(logbiomass=rep(1,times=2),
		comp=factor(rep(0:1,each=2)),herb=factor(rep(c("Ambient","Reduced"),2),levels=levels(BTseeds$herb)))

mm.slope = model.matrix(as.formula(BTseeds.models[[1]]),data=nd.slope)

mm.slope[,c(1,3:4,7)] = 0

int.results = array(NA,dim=c(length(BTseeds.fits),nrow(mm.intercept)))
int.se = matrix(NA,nrow=length(BTseeds.fits),ncol=nrow(mm.intercept))
slope.results = array(NA,dim=c(length(BTseeds.fits),nrow(mm.slope)))
slope.se = matrix(NA,nrow=length(BTseeds.fits),ncol=nrow(mm.intercept))

for (i in 1:length(BTseeds.fits)){
	fef = rep(0,times=length(coef.names))
	names(fef) = coef.names
	vcf = matrix(0,nrow=length(coef.names),ncol=length(coef.names))
	rownames(vcf) = coef.names
	colnames(vcf) = coef.names

	fe = fixef(BTseeds.fits[[i]])
	fef[names(fe)] = fe
	vcf[names(fe),names(fe)] = as.matrix(vcov(BTseeds.fits[[i]]))
	int.results[i,] = mm.intercept %*% fef
	slope.results[i,] = mm.slope %*% fef
	for (j in 1:nrow(mm.intercept)){
		pick = which(mm.intercept[j,]>0)
		int.se[i,j] = sqrt(sum(vcf[pick,pick]))
		pick = which(mm.slope[j,]>0)
		slope.se[i,j] = sqrt(sum(vcf[pick,pick])) # make it the se not the variance
	}
	
}

int.avg = apply(int.results,2,function(x)sum(x*wbic,na.rm=TRUE))
t1 = sweep(int.results,2,int.avg)
t2 = sqrt(int.se^2 + t1^2)
int.avg.se = apply(t2,2,function(x)sum(wbic*x,na.rm=TRUE))

slope.avg = apply(slope.results,2,function(x)sum(x*wbic,na.rm=TRUE))
t1 = sweep(slope.results,2,slope.avg)
t2 = sqrt(slope.se^2 + t1^2)
slope.avg.se = apply(t2,2,function(x)sum(wbic*x,na.rm=TRUE))


#comp = 0 refers to "Medium", and com = 1 refers to "High"

nd.match=nd.intercept
nd.match[,2] = c(rep("Medium",2),rep("High",2))
for (i in 1:4) {   IPM.parameters$seed.intercept[IPM.parameters$species=="BT" & IPM.parameters$comp==nd.match[i,"comp"] & IPM.parameters$herb==nd.match[i,"herb"]]=int.avg[i]}
for (i in 1:4) {IPM.parameters$seed.intercept.se[IPM.parameters$species=="BT" & IPM.parameters$comp==nd.match[i,"comp"] & IPM.parameters$herb==nd.match[i,"herb"]]=int.avg.se[i]}
for (i in 1:4) {       IPM.parameters$seed.slope[IPM.parameters$species=="BT" & IPM.parameters$comp==nd.match[i,"comp"] & IPM.parameters$herb==nd.match[i,"herb"]]=slope.avg[i]}
for (i in 1:4) {    IPM.parameters$seed.slope.se[IPM.parameters$species=="BT" & IPM.parameters$comp==nd.match[i,"comp"] & IPM.parameters$herb==nd.match[i,"herb"]]=slope.avg.se[i]}

write.csv(IPM.parameters,"IPM.Parameters.csv",row.names = FALSE)



#==================================================================
#Tall Thistle: data from Leland
#==================================================================

TTseeds2006=read.csv("Tall Thistle 2006 Updated Dec 2010.csv",header = T,na.strings=c("NA","na"))
TTseeds2007=read.csv("Tall Thistle 2007 Updated Dec 2010.csv",header = T,na.strings=c("NA","na"))
TTseeds=rbind(TTseeds2006,TTseeds2007)
TTseeds=TTseeds[-467,]

#converting variables in "bolt" data to numeric/categorical variables
TTseeds$Comp = factor(TTseeds$Comp,levels=levels(grow.tall$Comp))
TTseeds$Herb = factor(TTseeds$Herb,levels=levels(grow.tall$Herb))
TTseeds$site = as.factor(TTseeds$site)
TTseeds$year = as.factor(TTseeds$year)
TTseeds$site.year = paste(TTseeds$site,BTseeds$year,sep=".")
TTseeds$Topo=as.factor(TTseeds$Topo)

hist(TTseeds$total.seeds)
table(TTseeds$Comp,TTseeds$Herb)

TTseeds = TTseeds[complete.cases(TTseeds[,12]),]
plot(total.seeds~No.Lvs, data = TTseeds[TTseeds$year==2006,])
points(total.seeds~No.Lvs, data = TTseeds[TTseeds$year==2007,],col = "red")
TTseeds$logseeds = log(TTseeds$total.seeds+1)
TTseeds$logleaves = log(TTseeds$No.Lvs)
plot(logseeds~logleaves, data = TTseeds[TTseeds$year==2006,])
points(logseeds~logleaves, data = TTseeds[TTseeds$year==2007,],col = "red")

# now predict veg.bio07 for this dataset. 
# assumes that growth models have been fitted, check for appropriate models
if (!exists("tb.fits")) stop("fit growth models first!")

# now using model averaged predictions
nd = TTseeds
nd$no.leaves08  = nd$No.Lvs
nd$bolt08 = factor(1,levels=levels(grow.bull$bolt08))
results = matrix(NA,nrow=nrow(TTseeds),ncol=length(tb.fits))
for (i in 1:length(tb.fits)){
	# generate predictions from a lmer model using the data from inside the model, mostly. 
	fm = tb.fits[[i]]
	results[,i] = predict(fm,newdata=nd,re.form=~0)
}
# model averaged prediction for log.veg.bio parameter
# make sure we have the right weights vector!!
library(lme4)
bic = unlist(sapply(tb.fits,BIC)		)
dbic = bic-min(bic)
wbic = exp(-dbic/2)/sum(exp(-dbic/2))
TTseeds[,"log.veg.bio"] = apply(results,1,function(x)sum(x*wbic))

#-----------------------
#Probability of zero seed production under reduced herbivory is zero
TTseeds[,"seedproducedyesno"]= ifelse (TTseeds$total.seeds>0,1,0)

zeroseeds.mm0 = glmer(seedproducedyesno~log.veg.bio + Herb +  
		log.veg.bio:Herb + (1|site.year)+(1|Topo),  family=binomial,data=TTseeds)
zeroseeds.mm1 = glmer(seedproducedyesno~log.veg.bio + Herb +  
		log.veg.bio:Herb + (1|site.year),  family=binomial,data=TTseeds)
zeroseeds.mm2 = glmer(seedproducedyesno~log.veg.bio + Herb +  
		log.veg.bio:Herb + (1|Topo),  family=binomial,data=TTseeds)
zeroseeds.mm3 = glmer(seedproducedyesno~log.veg.bio + Herb +  
		log.veg.bio:Herb + (Topo|site.year),  family=binomial,data=TTseeds)
zeroseeds.mm4 = glmer(seedproducedyesno~log.veg.bio + Herb +  
		log.veg.bio:Herb + (1|site.year)+(Topo|site.year),  family=binomial,data=TTseeds)
zeroseeds.mm5 = glmer(seedproducedyesno~log.veg.bio + Herb +  
		log.veg.bio:Herb + (1|Topo)+(Topo|site.year),  family=binomial,data=TTseeds)

sapply(list(zeroseeds.mm0,zeroseeds.mm1,zeroseeds.mm2,zeroseeds.mm3,zeroseeds.mm4,zeroseeds.mm5),BIC)
# ok, (Topo|site.year) effect only
zeroseeds.models = list("~log.veg.bio + Herb + log.veg.bio:Herb",
				"~log.veg.bio + Herb",
				"~log.veg.bio",
				"~Herb",
				"~1")

zeroseeds.fits=lapply(zeroseeds.models,function(ff)glmer(paste("seedproducedyesno",ff,"+ (Topo|site.year)"),data=TTseeds,family=binomial))
bic = unlist(sapply(zeroseeds.fits,BIC)		)
dbic = bic-min(bic)
wbic = exp(-dbic/2)/sum(exp(-dbic/2))
data.frame(unlist(zeroseeds.models),dbic,wbic)[order(dbic),]

loglik=unlist(sapply(zeroseeds.fits,logLik))
k=sapply(zeroseeds.fits,function(mm)length(fixef(mm)))+2
stat.res=data.frame(unlist(zeroseeds.models),k,loglik,dbic,wbic)[order(dbic),]
write.csv(stat.res,"zeroseeds_res_TT.csv")


# make a plot
nd = data.frame(log.veg.bio=rep(seq(-0.1,5.1,0.1),times=2),
		Herb=factor(rep(c("Ambient","Reduced"),each=53),levels=levels(TTseeds$Herb)))

results = matrix(NA,nrow=nrow(nd),ncol=length(zeroseeds.fits))
for (i in 1:length(zeroseeds.fits)){
	# generate predictions from a glmer model using the data from inside the model, mostly. 
	fm = zeroseeds.fits[[i]]
	results[,i] = predict(fm,newdata=nd,re.form=~0)
}
# model averaged prediction for logbiomass parameter
nd[,"nonzero"] = apply(results,1,function(x)sum(x*wbic))

xyplot(nonzero~log.veg.bio|Herb,data=nd,type="l",
panel=function(...){
		panel.xyplot(...)
		pn = which.packet()
		#ltext(0,0.5,pn[1])
		pick = TTseeds$Herb == levels(TTseeds$Herb)[pn[1]] 			
		ticks = TTseeds[pick,c("log.veg.bio","viable.seeds")]
		panel.rug(ticks[ticks$viable.seeds=="N",1])
		panel.rug(ticks[ticks$viable.seeds=="Y",1],y=NULL, regular=FALSE)
	},
ylim=c(0,1))

coef.names = names(fixef(zeroseeds.fits[[1]]))
coef.results = matrix(NA,nrow=length(zeroseeds.fits),ncol=length(coef.names))
coef.se = matrix(NA,nrow=length(zeroseeds.fits),ncol=length(coef.names))
for (i in 1:length(zeroseeds.fits)){
	coef.results[i,] = sapply(coef.names,function(cf)fixef(zeroseeds.fits[[i]])[cf],USE.NAMES=FALSE)
	se = sqrt(diag(vcov(zeroseeds.fits[[i]])))
	names(se) = names(fixef(zeroseeds.fits[[i]]))
	coef.se[i,] = sapply(coef.names,function(cf)se[cf])
}

coef.avg = apply(coef.results,2,function(x)sum(x*wbic,na.rm=TRUE))
t1 = sweep(coef.results,2,coef.avg)
t2 = sqrt(coef.se^2 + t1^2)
coef.avg.se = apply(t2,2,function(x)sum(wbic*x,na.rm=TRUE))

nd.intercept = data.frame(log.veg.bio=rep(0,times=2),
		Herb=factor(c("Ambient","Reduced"),levels=levels(TTseeds$Herb)))
mm.intercept = model.matrix(as.formula(zeroseeds.models[[1]]),data=nd.intercept)

nd.slope = data.frame(log.veg.bio=rep(1,times=2),
		Herb=factor(c("Ambient","Reduced"),levels=levels(TTseeds$Herb)))

mm.slope = model.matrix(as.formula(zeroseeds.models[[1]]),data=nd.slope)

mm.slope[,c(1,3)] = 0

int.results = array(NA,dim=c(length(zeroseeds.fits),nrow(mm.intercept)))
int.se = matrix(NA,nrow=length(zeroseeds.fits),ncol=nrow(mm.intercept))
slope.results = array(NA,dim=c(length(zeroseeds.fits),nrow(mm.slope)))
slope.se = matrix(NA,nrow=length(zeroseeds.fits),ncol=nrow(mm.intercept))

for (i in 1:length(zeroseeds.fits)){
	fef = rep(0,times=length(coef.names))
	names(fef) = coef.names
	vcf = matrix(0,nrow=length(coef.names),ncol=length(coef.names))
	rownames(vcf) = coef.names
	colnames(vcf) = coef.names

	fe = fixef(zeroseeds.fits[[i]])
	fef[names(fe)] = fe
	vcf[names(fe),names(fe)] = as.matrix(vcov(zeroseeds.fits[[i]]))
	int.results[i,] = mm.intercept %*% fef
	slope.results[i,] = mm.slope %*% fef
	for (j in 1:nrow(mm.intercept)){
		pick = which(mm.intercept[j,]>0)
		int.se[i,j] = sqrt(sum(vcf[pick,pick]))
		pick = which(mm.slope[j,]>0)
		slope.se[i,j] = sqrt(sum(vcf[pick,pick])) # make it the se not the variance
	}
	
}

int.avg = apply(int.results,2,function(x)sum(x*wbic,na.rm=TRUE))
t1 = sweep(int.results,2,int.avg)
t2 = sqrt(int.se^2 + t1^2)
int.avg.se = apply(t2,2,function(x)sum(wbic*x,na.rm=TRUE))

slope.avg = apply(slope.results,2,function(x)sum(x*wbic,na.rm=TRUE))
t1 = sweep(slope.results,2,slope.avg)
t2 = sqrt(slope.se^2 + t1^2)
slope.avg.se = apply(t2,2,function(x)sum(wbic*x,na.rm=TRUE))


# We only have data on Comp = "High"
IPM.parameters$zero.seed.intercept[IPM.parameters$species=="TT" & IPM.parameters$comp=="High" & IPM.parameters$herb=="Ambient"]=int.avg[1]
IPM.parameters$zero.seed.intercept[IPM.parameters$species=="TT" & IPM.parameters$comp=="High" & IPM.parameters$herb=="Reduced"]=int.avg[2]
IPM.parameters$zero.seed.intercept.se[IPM.parameters$species=="TT" & IPM.parameters$comp=="High" & IPM.parameters$herb=="Ambient"]=int.avg.se[1]
IPM.parameters$zero.seed.intercept.se[IPM.parameters$species=="TT" & IPM.parameters$comp=="High" & IPM.parameters$herb=="Reduced"]=int.avg.se[2]

IPM.parameters$zero.seed.slope[IPM.parameters$species=="TT" & IPM.parameters$comp=="High" & IPM.parameters$herb=="Ambient"]=slope.avg[1]
IPM.parameters$zero.seed.slope[IPM.parameters$species=="TT" & IPM.parameters$comp=="High" & IPM.parameters$herb=="Reduced"]=slope.avg[2]
IPM.parameters$zero.seed.slope.se[IPM.parameters$species=="TT" & IPM.parameters$comp=="High" & IPM.parameters$herb=="Ambient"]=slope.avg.se[1]
IPM.parameters$zero.seed.slope.se[IPM.parameters$species=="TT" & IPM.parameters$comp=="High" & IPM.parameters$herb=="Reduced"]=slope.avg.se[2]



#===============================
#Remove zero seed production
remove=which(TTseeds$total.seeds==0)

TTseeds.mm0=lmer(logseeds~log.veg.bio + herb + log.veg.bio:herb +
		(1|site.year)+(1|Topo),  data=TTseeds[-remove,])
TTseeds.mm1=lmer(logseeds~log.veg.bio + herb + log.veg.bio:herb +
		(1|site.year),  data=TTseeds[-remove,])
TTseeds.mm2=lmer(logseeds~log.veg.bio + herb + log.veg.bio:herb +
		(1|Topo),  data=TTseeds[-remove,])
TTseeds.mm3=lmer(logseeds~log.veg.bio + herb + log.veg.bio:herb +
		(1|site.year)+(Topo|site.year),  data=TTseeds[-remove,])
TTseeds.mm4=lmer(logseeds~log.veg.bio + herb + log.veg.bio:herb +
		(1|site.year)+(Topo|site.year),  data=TTseeds[-remove,])
TTseeds.mm5=lmer(logseeds~log.veg.bio + herb + log.veg.bio:herb +
		(1|Topo)+(Topo|site.year),  data=TTseeds[-remove,])

sapply(list(TTseeds.mm0,TTseeds.mm1,TTseeds.mm2,TTseeds.mm3,TTseeds.mm4,TTseeds.mm5),AIC)
Res1=resid(TTseeds.mm2)
plot(fitted(TTseeds.mm2),Res1)

qqnorm(Res1)
qqline(Res1)

TTseeds.models = list("~log.veg.bio + Herb  + log.veg.bio:Herb",
		"~log.veg.bio + Herb",
		"~Herb",
		"~log.veg.bio",
		"~1")
TTseeds.fits = lapply(TTseeds.models,function(ff)lmer(paste("logseeds",ff,"+ (1|Topo)"),data=TTseeds[-remove,]))

bic = unlist(sapply(TTseeds.fits,BIC)		)
dbic = bic-min(bic)
wbic = exp(-dbic/2)/sum(exp(-dbic/2))
data.frame(unlist(TTseeds.models),dbic,wbic)[order(dbic),]
loglik=unlist(sapply(TTseeds.fits,logLik))
k=sapply(TTseeds.fits,function(mm)length(fixef(mm)))+2
stat.res=data.frame(unlist(TTseeds.models),k,loglik,dbic,wbic)[order(dbic),]
write.csv(stat.res,"Seeds_res_TT.csv")

# make a plot
nd = data.frame(log.veg.bio=rep(seq(0.1,5.1,0.1),times=2),
		Herb=factor(rep(c("Ambient","Reduced"),each=51),levels=levels(TTseeds$Herb)))


results = matrix(NA,nrow=nrow(nd),ncol=length(TTseeds.fits))
for (i in 1:length(TTseeds.fits)){
	# generate predictions from a glmer model using the data from inside the model, mostly. 
	fm = TTseeds.fits[[i]]
	results[,i] = predict(fm,newdata=nd,re.form=~0)
}
# model averaged prediction for veg.bio07 parameter

nd[,"noseeds"] = apply(results,1,function(x)sum(x*wbic))

xyplot(noseeds~log.veg.bio|Herb,data=nd,type="l",
panel=function(...){
		panel.xyplot(...)
		pn = which.packet()
		pick = TTseeds$Herb == levels(TTseeds$Herb)[pn[1]] 
		ticks = TTseeds[pick,c("log.veg.bio","logseeds")]
		lpoints(ticks$log.veg.bio,ticks$logseeds)
	}
)


#------------------
coef.names = names(fixef(TTseeds.fits[[1]]))
coef.results = matrix(NA,nrow=length(TTseeds.fits),ncol=length(coef.names))
coef.se = matrix(NA,nrow=length(TTseeds.fits),ncol=length(coef.names))
for (i in 1:length(TTseeds.fits)){
	coef.results[i,] = sapply(coef.names,function(cf)fixef(TTseeds.fits[[i]])[cf],USE.NAMES=FALSE)
	se = sqrt(diag(vcov(TTseeds.fits[[i]])))
	names(se) = names(fixef(TTseeds.fits[[i]]))
	coef.se[i,] = sapply(coef.names,function(cf)se[cf])
}

coef.avg = apply(coef.results,2,function(x)sum(x*wbic,na.rm=TRUE))
t1 = sweep(coef.results,2,coef.avg)
t2 = sqrt(coef.se^2 + t1^2)
coef.avg.se = apply(t2,2,function(x)sum(wbic*x,na.rm=TRUE))

nd.intercept = data.frame(log.veg.bio=rep(0,times=2),
		   Herb=factor(c("Ambient","Reduced"),levels=levels(TTseeds$Herb)))

mm.intercept = model.matrix(as.formula(TTseeds.models[[1]]),data=nd.intercept)

nd.slope = data.frame(log.veg.bio=rep(1,times=2),
		Herb=factor(c("Ambient","Reduced"),levels=levels(TTseeds$Herb)))

mm.slope = model.matrix(as.formula(TTseeds.models[[1]]),data=nd.slope)

mm.slope[,c(1,3)] = 0

int.results = array(NA,dim=c(length(TTseeds.fits),nrow(mm.intercept)))
int.se = matrix(NA,nrow=length(TTseeds.fits),ncol=nrow(mm.intercept))
slope.results = array(NA,dim=c(length(TTseeds.fits),nrow(mm.slope)))
slope.se = matrix(NA,nrow=length(TTseeds.fits),ncol=nrow(mm.intercept))

for (i in 1:length(TTseeds.fits)){
	fef = rep(0,times=length(coef.names))
	names(fef) = coef.names
	vcf = matrix(0,nrow=length(coef.names),ncol=length(coef.names))
	rownames(vcf) = coef.names
	colnames(vcf) = coef.names

	fe = fixef(TTseeds.fits[[i]])
	fef[names(fe)] = fe
	vcf[names(fe),names(fe)] = as.matrix(vcov(TTseeds.fits[[i]]))
	int.results[i,] = mm.intercept %*% fef
	slope.results[i,] = mm.slope %*% fef
	for (j in 1:nrow(mm.intercept)){
		pick = which(mm.intercept[j,]>0)
		int.se[i,j] = sqrt(sum(vcf[pick,pick]))
		pick = which(mm.slope[j,]>0)
		slope.se[i,j] = sqrt(sum(vcf[pick,pick])) # make it the se not the variance
	}
	
}

int.avg = apply(int.results,2,function(x)sum(x*wbic,na.rm=TRUE))
t1 = sweep(int.results,2,int.avg)
t2 = sqrt(int.se^2 + t1^2)
int.avg.se = apply(t2,2,function(x)sum(wbic*x,na.rm=TRUE))

slope.avg = apply(slope.results,2,function(x)sum(x*wbic,na.rm=TRUE))
t1 = sweep(slope.results,2,slope.avg)
t2 = sqrt(slope.se^2 + t1^2)
slope.avg.se = apply(t2,2,function(x)sum(wbic*x,na.rm=TRUE))


# We only have data on Comp = "High"
IPM.parameters$seed.intercept[IPM.parameters$species=="TT" & IPM.parameters$comp=="High" & IPM.parameters$herb=="Ambient"]=int.avg[1]
IPM.parameters$seed.intercept[IPM.parameters$species=="TT" & IPM.parameters$comp=="High" & IPM.parameters$herb=="Reduced"]=int.avg[2]
IPM.parameters$seed.intercept.se[IPM.parameters$species=="TT" & IPM.parameters$comp=="High" & IPM.parameters$herb=="Ambient"]=int.avg.se[1]
IPM.parameters$seed.intercept.se[IPM.parameters$species=="TT" & IPM.parameters$comp=="High" & IPM.parameters$herb=="Reduced"]=int.avg.se[2]

IPM.parameters$seed.slope[IPM.parameters$species=="TT" & IPM.parameters$comp=="High" & IPM.parameters$herb=="Ambient"]=slope.avg[1]
IPM.parameters$seed.slope[IPM.parameters$species=="TT" & IPM.parameters$comp=="High" & IPM.parameters$herb=="Reduced"]=slope.avg[2]
IPM.parameters$seed.slope.se[IPM.parameters$species=="TT" & IPM.parameters$comp=="High" & IPM.parameters$herb=="Ambient"]=slope.avg.se[1]
IPM.parameters$seed.slope.se[IPM.parameters$species=="TT" & IPM.parameters$comp=="High" & IPM.parameters$herb=="Reduced"]=slope.avg.se[2]


write.csv(IPM.parameters,"IPM.Parameters.csv",row.names = FALSE)